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ABSTRACT 

Context. Recent work identified a growth barrier for dust coagulation that originates in the electric repulsion between 
colliding particles. Depending on its charge state, dust material may have the potential to control key processes towards 
planet formation such as MHD (magnetohydrodynamic) turbulence and grain growth which are coupled in a two-way 
process. 

Aims. We quantify the grain charging at different stages of disc evolution and differentiate between two very extreme 
cases: compact spherical grains and aggregates with fractal dimension Df = 2. 

Methods. Applying a simple chemical network that accounts for collisional charging of grains, we provide a semi- 
analytical solution. This allowed us to calculate the equilibrium population of grain charges and the ionisation fraction 
efficiently. The grain charging was evaluated for different dynamical environments ranging from static to non- stationary 
disc configurations. 

Results. The results show that the adsorption/desorption of neutral gas- phase heavy metals, such as magnesium, effects 
the charging state of grains. The greater the difference between the thermal velocities of the metal and the dominant 
molecular ion, the greater the change in the mean grain charge. Agglomerates have more negative excess charge on 
average than compact spherical particles of the same mass. The rise in the mean grain charge is proportional to iV 1 / 6 
in the ion-dust limit. We find that grain charging in a non-stationary disc environment is expected to lead to similar 
results. 

Conclusions. The results indicate that the dust growth and settling in regions where the dust growth is limited by the 
so-called "electro-static barrier" do not prevent the dust material from remaining the dominant charge carrier. 
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1. Introduction 

Small dust particles are regarded as the key ingredient 
to control the planet formation process in protoplanetary 
discs. At early stages of planet formation, submicron-sized 
compact particles are thought to grow towards larger but 
fluffy agglomerates. At mm to cm sizes these agglomerates 
are supposed to be compacted, whereby they loose their 
fractal structure. Because of the rise in the mass-to-surface 
ratio associated with the compaction, agglomerates tend to 
be more concentrated locally. However, there are a number 
of potential processes that may control grain growth, one 
of which is disc turbulence. The magnetorotational insta- 
bility - MRI - (Balbus & Hawley 1991; Hawley & Balbus 
1991) has been shown to be robust in generating turbu- 
lence in Keplerian discs. MHD turbulence provides the in- 
ternal stress required for mass accretion, the rate of which 
is constrained by observations. Observations of young stars 
indicate that discs usually show signatures of active gas ac- 
cretion onto the central star with mass flow rates of about 
10- 8±1 Mgyr -1 (e.g., Sicilia-Aguilar et al. 2004). Recent 
observations from discs around young stars also set con- 
straints on the turbulent linewidth and provide support for 
subsonic turbulence in discs (Hughes et al. 2011). However, 
numerical studies (e.g. Fleming & Stone 2003) on MRI- 
driven MHD turbulence have identified locations within the 
planet forming regions, the so-called "dead zones", where 



the coupling between the gas and the magnetic field is not 
sufficient to maintain MHD turbulence. Fleming & Stone 
also showed that a low Reynolds stress can be maintained 
in the dead zone, such that low levels of accretion are sus- 
tained there. That is why dead zones have been considered 
to advance the planet formation process. Dead zones pro- 
vide a safe environment for grain growth and for planetes- 
imals (Gressel et al. 2011). However, the presence of small 
grains in the weakly ionized dead zones leads to significant 
changes in the abundances of the charge carriers in the gas 
phase and therefore affects the MRI. In other words, MHD 
turbulence and grain growth are coupled in a two-way pro- 
cess. 

Grain charging may also effect the dust growth: 
Okuzumi (2009) pointed out that the so-called "electro- 
static barrier" may inhibit dust coagulation in planet form- 
ing regions. In a subsequent study, Okuzumi et al. (2011a) 
extended the analysis including the dust size distribution. 
They found that under certain conditions the dust under- 
goes bimodal growth. While the small aggregates sweep 
up free electrons, the large aggregates continue to grow. 
However, the growing stalls if the electrostatic repulsion be- 
comes strong and the collisions are driven by Brownian mo- 
tion. The results of Okuzumi et al. (2011b) indicate that un- 
der minimum mass solar nebula conditions the dust growth 
halts at aggregate sizes beyond [4-10 -5 , 3-10 -2 ] cm depend- 
ing on the radial position. Conventional chemical models 
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studying the dust-grain chemistry in protoplanetary discs 
account for up to \Z\ = 3 grain charges (e.g. Sano et al. 
2000). More recently, the range of grain charges was ex- 
tended considerably for various purposes. Focusing on sur- 
face layers of T-Tauri and transitional discs, Perez-Becker & 
Chiang (2011) examined the charge distribution on grains 
covering grain charges [Z m i„, 2 m „] = [—200, 200]. Okuzumi 
(2009) calculated the charge distribution on dust aggre- 
gates. He demonstrated that aggregates made of N = 10 10 
monomers can carry excess charges of about Z m i n ~ — 10 5 . 
To account for higher grain charges we therefore improved 
our simple chemical network introduced in Ilgner & Nelson 
(2006a) allowing [Z m i n , Z max ] to become free of choice and 
appropriate to the conditions applied. 

In this paper we examine the grain charging for var- 
ious stages of disc evolution ranging from static to non- 
stationary disc environment. Linking the grain charging 
with the gas-dust dynamics provides a natural environment 
to mimic the variation of the dust-to-gas ratio Ed/S g . The 
grain charging is assumed to originate in collisional charg- 
ing processes between grains, electrons, and gas-phase ions. 
We include X-ray ionisation from the central star as the 
primary source for ionisation and consider, to some extent, 
the contributions from cosmic ray particles. We further- 
more assume an equilibrium distribution of dust material 
in vertical direction balancing turbulent stirring and sedi- 
mentation. In particular, we apply the model of Takeuchi 
& Lin (2005) to simulate the dynamics of the gas-dust disc. 
One of our primary goals is to compare the grain charging 
obtained for compact spheres and dust agglomerates and 
the mean grain charge in particular. 

We have investigated the effect that thermal adsorp- 
tion/desorption of metals has on grain charging using two 
different chemical networks. In comparison with the val- 
ues obtained by switching off the thermal adsorption of 
metals, the modified Oppenheimcr-Dalgarno model pro- 
duces lower, the Umebayashi-Nakano model higher values 
for the mean grain charge. This originates in differences 
in the thermal velocities of the dominant gas-phase ion: 
v HGO+ < v Mg+ < %+• I n une with expectations, we find 
that dust agglomerates have a higher charge-to-mass ra- 
tio carrying more charges than the corresponding compact 
spherical particles. Considering stationary disc configura- 
tions, we observe that grain charging on dust agglomerates 
is always associated with the so-called "ion-dust regime", 
where the charge balance is mainly maintained by nega- 
tively charged grains and positively charged gas-phase ions. 
(The ion-electron regime is characterised by the balance be- 
tween gas-phase ions and electrons.) In particular, we show 
that the fractional abundances of the charge carriers are 
independent of the number N of constituent monomers. 
Concerning compact spherical dust particles, we identified 
regions z/h g > associated with the ion-electron regime. 
Another conclusion of our work is that - for the parame- 
ter range considered (i.e., compact radii a* < 10~ 3 cm and 
agglomerates with characteristic radii a c < 10~ 2 cm with 
N < 10 6 , respectiveljQ) - grain charging in a non-stationary 
disc environment is expected to lead to similar results. We 
also present semi-analytical solutions for both agglomerates 
and compact spheres that allow us to determine the equi- 
librium distribution of <Z>, V <AZ 2 >, and the fractional 



1 The characteristic and compact radius is defined in Eqs. (|16p 
and (|18[l . respectively. 



abundances of the charge carriers a priori. 

The paper is organised as follows. In Sect. 2 we discuss 
the disc model we apply. In Sect. 3 we introduce the chem- 
ical network that we used to compare our results with the 
reaction network Okuzumi (2009) applied. Key problems 
concerning the ionisation rates are discussed in Sect. 4 while 
the numerical methods applied are described in Sect. 5. In 
Sect. 6 we present the results of our model, and finally in 
Sect. 7 we summarise the key findings of our study. 



2. Disc model 

The geometry of the disc model is described in cylindrical 
coordinates (R,(p,z), where R, tp, and z denote the radial 
distance to the point of interest, the azimuthal angle, and 
the height. Here, we simply considered a two-dimensional 
(2D) geometry by assuming axial symmetry, i.e., d/dip = 0. 

Treating the disc as a two-component fluid of dust par- 
ticles and gas, we assumed that the circumstellar gas disc 
is turbulent and accreting onto the central star. We also 
assumed that the gas is primarily heated by stellar radia- 
tion rather than viscous dissipation. We supposed that the 
mean gas motion is characterised by a subsonic flow in R, 
supersonic flow in ip, while the gas density is characterised 
by hydrostatic equilibrium in z-direction. 

The underlying disc model we applied is based on the 
assumption of a two-component fluid that allows the dust 
particles and the gas to follow different velocity fields 
and u g , respectively. In particular, we considered dust par- 
ticles for which the two-component fluid is best described 
in the so-called "free molecular approximation", e.g., ao <C 
1.44 i?ai/ 4 cm under minimum mass solar nebula condition, 
(ao denotes the dust radius while i? au is the orbital radius in 
units of AU.) Since the mean free path of the gas molecules 
exceeds the dust particle sizes, the collisions between the 
gas and the dust particles are much more frequent than the 
collisions between gas particles. We moreover considered 
disc stages for which the gas density g g dominates over the 
density ga of the dust. The gas is therefore assumed to be 
decoupled from the dust disc. However, the gas drag may al- 
ter the velocity Vd of the dust particles. Assuming iij < c s , 
the drag coefficient j3 is characterised by (Epstein, 1924) 



P = gather , (1) 

where Uth denotes the thermal velocity. We can addition- 
ally specify the friction or stopping time t$ — m//3, which 
will turn out to be an important diagnostic quantity, m and 
<j refer to the mass of the dust particle and its projected 
surface area, respectively. The friction time basically char- 
acterises the time scale needed to couple the dust dynamics 
with the gas, which is why ts is also called stopping time. 
Instead of is, we used the dimensionless stopping time Ts 
defined as the ratio between is and the dynamical time 
scale tip. 

We studied the grain charging for various stages of disc 
evolution ranging from static to non-stationary disc envi- 
ronment. A description of the gas-dust dynamics applied is 
given in Subsec. I2TT1 while Subsec. |2"T21 summarises the local 
disc structure. 
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2.1. Dynamics of the gas-dust disc 

The extreme case of disc dynamics is associated with the 
so-called "minimum mass solar nebula" (MMSN) model, 
which was proposed by Hayashi (1981). He adopted the 
following radial profiles of the dust and gas surface densities 
E d and E K for R < 36 AU 



E g = 1.7 • 10 3 i?- 3 / 2 gcm- 



E d 



'7.1 i?au /2 gcm- 2 
30 i?au 3/2 gcm- 2 



if 0.35 < i? au < 2.7 
if 2.70 < i? au < 36 . 



(2) 
(3) 



We recall that Hayashi derived the power law above assum- 
ing no substantial transport of dust material. Although re- 
cent simulations of the radial migration of planets indicate 
that the radial profile of E g (as well as the gas tempera- 
ture profile) differs significantly from the MMSN model (cf 
Desch 2007), we opted for Hayashi's MMSN model. That 
is because the assumed power index for the MMSN model 
nicely corresponds to a static solution of the disc model of 
Takeuchi & Lin (2002), see below for more information. The 
grain charging obtained for that particular environment is 
discussed in Sec. |3] 

Takeuchi & Lin (2005) extended their disc model to 
study the dynamical evolution of the gas-dust disc. The 
corresponding set of continuity equations is 



SE g 1 d ,__ _ . 
-W + RdR^^ 




0, 



(4) 
(5) 



where VR, g and T)R t d are the vertically integrated velocities 
of the gas and the dust particles. The (ordinary) diffusion 
is described by the phenomenological equation, which is 
well-known as Fick's first law for a binary system. The cor- 
responding mass flow is approximated by 



F dis 



v_d_ 

•S^dR 



Ed 

E„- 



(6) 



with the Schmidt number Sc = vjT> specifying the ratio of 
the rate of angular momentum transport to the rate of tur- 
bulent transport of grain particles. We used the a prescrip- 
tion for the viscous stress, such that the kinematic viscosity 
is v = ac s h g . 

The interpretation of the integrated velocities ^_R, g and 
w_R, d is crucial for our understanding of the species' trans- 
port in our particular model. VR ig and VR.d refer rather to 
the velocities associated with the transport of the surface 
densities E g and E d than to the velocities the individual 
species are transported with. For the gas disc, the velocity 
VR, g can be obtained by combining the equation of angular 
momentum and the continuity equation 



v R,s 



RV 2 E e dR 



(7) 



Inserting equation Q into the continuity equation for the 
gas surface density, we obtain the well-known diffusion 
equation reviewed, e.g., in Pringle (1981). 

Assuming power laws for the pressure scale height h g , 
the sound speed c s , and the local gas density g g (see Subsec. 
2.21 below). Takeuchi & Lin (2002) have presented steady- 
state solutions of Eqs. (0]) - ©. Evaluating the mass fluxes 



.REgt^g and RE^vr^, they obtained E g cx R 3 / 2 for 
VR tg /c s — and 



Eg — E i? au 

E d = M d / (2ttRv ra ) 



(8) 
(9) 



for VR tg /c s where the mass flux Md of dust material 
enters as a free parameter, vr d 



VR,d = / QdVR,ddz 

^d 



(10) 



can an be calculated semi-analytically (see Takeuchi & Lin, 
2002) assuming 



S v R.g - ^K.mid 

Ts+T- 1 



(11) 



Here, r\ and VK.mid refer to the ratio of the pressure gradi- 
ent to the gravity and Keplerian velocity at disc midplane, 
respectively. 

However, apart from steady-state conditions it is diffi- 
cult to find a reliable approximation for vr^&. We followed 
the instruction of Takeuchi & Lin (2005) to calculate vr^ 
but we are aware of the caveats inherent in their approach. 
This problem arises because a general analytical expression 
for the corresponding radial component VR tg of the gas ve- 
locity field does not exist. In order to tackle this problem 
we simply approximate Vra by 



Vra 



TFK.mid 



evaluating Tg and 77 at disc midplane: 



V 



R 

3 7T TO 

Ts = 8E^7 



<91nE„ 

& 

dlnR 



(12) 



(13) 



(14) 



where q denotes the power index associated with the gas 
pressure scale height h g (see next subsection). 

As we will demonstrate below, the dynamics of the gas- 
dust disc is controlled by the stopping time Tg. Because 
of Ts oc m/cr the stopping time relies on the particular 
topology of the dust particles considered. We studied two 
extreme cases of grain topology: compact spherical grains 
and irregular, very porous dust agglomerates. We begin 
addressing that particular question for compact spherical 
grains before discussing the dust agglomerates and its ef- 
fect on the gas-dust dynamics. 

Case A : Compact spherical dust grains 

Considering compact spherical grain particles of material 
bulk density g p and a radius do, the stopping time at disc 
midplane z/h g = is given by 



Ts = 



2E ff ' 



Takeuchi & Lin (2005) studied in detail the effect of com- 
pact spherical grains on the dynamics of the gas-dust disc. 
We took their fiducial model, which nicely demonstrates 
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Fig. 1. Radial profile of the surface densities of the gas and the dust at different time snaps t = 10 5 , 3 • 10 5 and 10 6 yr. 
The solid lines correspond to the dust surface density while the gas surface density is shown using the dashed lines. The 
(red-coloured) dotted line refers to S g - R- 1 and E d = M d / ! (2itR<v ra >) at t/t K = 0. \M d \ = KT lo M /yr, a = 10~ 3 , 
and g p = 1.0 gem -3 . The profiles shown in the left panel are obtained under the assumption of spherical grains with 
ao = 10~ 3 cm. The corresponding profiles obtained for agglomerates with -Dbcca = 2 and N = 10 6 monomers of 
do = 10 -5 cm are shown in the right panel. 



the inward migration of larger (spherical) dust particles to 
test our implementation of the numerical schemes applied. 
We also adopted boundary and initial conditions suitable 
for our particular modelJ3 To summarise, we set the inner 
boundary at i?i rm cr = 0.1 AU and replaced the zero torque 
boundary conditions with predefined analytical values for 
E g (i). In particular, we applied the self-similarity solution 

£ g = SqX" 1 exp{-aX/r} . (15) 

Lynden-Bell & Pringle (1974) obtained under the assump- 
tion i>t oc R with dimensionless variables X = R/Ro and 
r = t/to + 1. In the limit of X <C 1 the expression re- 
duces to S g oc A -1 , which represents the steady-state 
solution Eg oc discussed above. The most appropri- 
ate value for the constant a was applied to approximate 
£ g (t) in the outer disc regions. We replaced the initial as- 
sumption Ed/Sg = const by \Md\ = const. For the sake 
of convenience we applied the value of Takeuchi & Lin 
(2002) M d = 10 -10 M©yr -1 , which is one order of mag- 
nitude smaller than the corresponding mass flux of the gas 
M g = 37w£ g - 2.6 x lO^Mgyr" 1 . 

Because of Tg oc ao, the dynamics of the dust disc de- 
pends on the size of the spherical grain particle considered. 
This may also have important consequences for the dust 
growth. The simplest approach to mimic the dust growth is 
to assume that the entire available dust mass is represented 
by monodisperse compact spherical particles of a given size 
a* . By changing a» one may trace different stages of dust 
growth. However, dust growth is supposed to be processed 
by coagulation leading to irregular dust agglomerates. The 
effect of dust agglomerates on the dust dynamics is dis- 
cussed below. 

Case B : Irregular porous dust agglomerates 

The grain particles are now regarded as dust agglomerates. 

2 We investigated these effects; the results obtained are dis- 
cussed in the appendix. 



In particular, we considered agglomerates associated with 
the so-called "ballistic cluster-cluster aggregation" (BCCA) 
model. The BCCA model describes the growth of aggre- 
gate^ through collisions of equal-sized aggregates as car- 
tooned in Fig. [3J We are aware of constraints associated 
with the BCCA growth process (see discussion in Okuzumi 
2009). However, for the purpose of our paper we simply as- 
sumed that BCCA agglomerates are formed independently 
of the disc environment applied. That is because the BCCA 
agglomerates are used as test particles to trace the grain 
charging for different disc environments. 

Simulations of dust growth have shown that for the 
BCCA model the characteristic radiufQ a c obeys a power 
law (e.g., Okuzumi et al. 2009) while the mean value for 
the projected surface area is well approximated by a linear 
function of N (cf Minato et al. 2006): 

a c = a N 1/DBCCA (16) 
ctbcca = (0.351A + 0.5667V°- 862 )7ra^ for N > 16. (17) 

N denotes the number of the constituent monomers forming 
the agglomerate (ao is the radius of each monomer). £>bcca 
is the fractal dimension of BCCA associated with a c and 
■Dbcca ~ 1.9 (Meakin 1991). We considered -Dbcca = 2. 
The projected surface area is ctbcca ~ ttOc = Nitclq, which 
is in line with the linear approximation of Eq. (|17p . The 
expression of the stopping time (fT4"]) then is 

7T g p a 



3 We simplified the terminology used in this paper. Because 
fractal agglomerates with fractal dimension D ~ 2 can be mod- 
elled applying the BCCA, agglomerates with D ~ 2 are referred 
to as "BCCA agglomerates". 

4 The characteristic radius a c serves to measure the size of 
the agglomerate. In particular, its weighting is different from 
the rms-definition of the gyration radius r s , cf. Okuzumi et al. 
(2009). 
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scaling exponent q = —1/2: 



Ballistic cluster-cluster aggregation (BCCA) 



o - -» * 




-icrcmers 



Fig. 2. Scheme representing the dust growth characterised 
as ballistic cluster-cluster aggregation. The agglomerates 
collide with a clone of each as marked with filled circles. 
The figure is adopted from Okuzumi et al. (2009). 



tig - «0,grt A U 

c s = co-R^u 

so that the gas density is defined by 
1 E„ ( 1 z 2 ~ 



/2tt K 



■ exp 



2 hi 



(19) 
(20) 

(21) 



Regarding the dust disc, we are primarily interested in con- 
ditions for which the mass flows associated with the set- 
tling and (vertical) turbulent mixing of the dust particles 
approach or attain a state of equilibrium. Taking the fast 
settling limit T s /a 3> (h g /R) 2 into account, the dust den- 
sity profile becomes 



Qd = 



1 Z ^s.mid 

Po.dexp < --— - be- 



2 ^ 



a 



Qddz, 



(22) 
(23) 



We note that for that particular model the stopping time 
is independent of the number of monomers considered. We 
conclude that the agglomerates of different sizes (i.e., differ- 
ent N) experience the same dynamical evolution contrast- 
ing the properties derived for compact spherical grains. 
We can relate the single agglomerate made of N monomers 
of size clq to a compact sphere of the same mass by 



cf Takeuchi & Lin (2002). Approximating the settling time 
scale by t sc< \ — (Tgilx) -1 , we confirm in the fast settling 
limit the validity of 



t z <C i S ed <C t v 



1 



1 



1 



(h g /R) 2 K ' 



(24) 



We considered that the gas is transparent to the emitted 
visible radiation of the protostar and neglected the contri- 
butions to the gas temperature from viscous dissipation. 
The temperature T g of the gas is then given by 



a* = agN 



1/3 



(18) 



which turns out to become an important diagnostic quan- 
tity for our calculations. According to Eq. (TTS|) . a compact 
sphere of a* = 10~ 3 cm corresponds to an agglomerate of 
N = 10 6 monomers with ao = 10 -5 cm. The surface den- 
sity profiles obtained for this dust agglomerate are shown in 
the right panel of Fig. [T] We observe that the dust is tightly 
coupled to the gas dynamics that expands because of tur- 
bulent mixing towards larger orbital radii R. The advective 
inward drift of the dust (which is the dominant transport 
process for compact spherical grains ao > 10 -2 cm) cannot 
compete with the diffusive mixing along R. 



2.2. Local disc structure 



R 



AU 



(25) 



with T = 280 K (Hayashi 1981). Regarding the temper- 
ature of the dust particles, we drastically simplified the 
physics by setting Td = T g , which is in line with pre- 
vious studies (Sano et al. 2000, Bai & Goodman 2009). 
However, we are aware of the effect stellar radiation has 
on dust particles at high altitude (i.e., the so-called "super- 
heated layer" ) and its consequences for the interior gas tem- 
perature as demonstrated by Chiang & Goldreich (1997). 
Chiang & Goldreich report that the superheated dust layer 
may change with R; for their particular model of the min- 
imum mass solar nebula, they showed that the location of 
the superheated layer changes from z / h g = 5 at R = 3 AU 
to z/h s = 4 at R = 100 AU. For altitudes z/h g < 3 we 
considered, that the dust particles may be located some 
distance below the superheated layer. 



We applied the conventional ansatz of the thin disc ap- 
proximation h s /R -C 1. We assumed that the gas disc is 
isothermal in z-direction. Then, we know that on a time 
scale t z — h g /c s — Q^ 1 the hydrostatic equilibrium is es- 
tablished while changes in the gas surface density E g are de- 
termined by the viscous time scale i„ = R 2 jv with t z <C t„. 
Hence the local structure of the gas disc is assumed to set 
up instantaneously during the viscous evolution. 

Following Takeuchi & Lin (2002), we assumed that the 
gas pressure scale height h g and the isothermal sound speed 
c s are simply defined by power laws associated with the 



3. Chemical model 

We briefly recall the constituents of a chemical model: its 
components (or species), the chemical reactions that cause 
time-dependent changes in the abundances of species, and 
the underlying kinetic equation for each component. 

Regarding the class of chemical species applied, we con- 
sidered gas-phase species, grain species, and mantle species. 
The latter are species adsorbed onto the grain surface and 
associated with the counterparts of the neutral gas-phase 
species. Using this nomenclature, we differentiate between 
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Table 1. List of reactions considered for the modified 
Oppcnheimcr-Dalgarno network. M and m denote the neu- 
tral metal and molecule component, represented by magne- 
sium and molecular hydrogen, while X = {Mg, H2}. kr, and 
&6 represent the chemical processes that are not covered by 
the Umebayashi-Nakano network. 



m +M + 
111 

M + hv 
m + + e~ 
gX 

X + grains 



fcl 


m+ 


+ 


M 




m+ 


+ 


e~ 


k. 


M+ 


+ 


e~ 


A'4 


m 






kr, 


X 


+ 


grains 


k(i 


gX 






k- 


x+ 


+ 


gr Z 



gX 

x 



.z+i 
.z+i 



for Z > 
else 



reactions that involve mantle species and those that do not. 
As we did in Ilgner & Nelson (2006a), we refer to the for- 
mer as mantle chemistry and the latter as grain chemistry. 
We considered two different types for the topology of grain 
particles: 

Type 1: The grain particles are regarded as dust agglom- 
erates made by dust growth. For our particular model, we 
consider agglomerates associated with the BCCA model. 
The individual agglomerates are also assigned to different 



grain charges with j 
Sec. ~ 



•Zmax- Details are given in 



Type 2: The grain particles are now characterised as spher- 
ical dust grains of a radius clq. The individual grain par- 
ticles are assigned to different grain charges with j = 
Z m i n , . . . , -Zmax. These models are discussed in Sec. 13.31 

For the chemical reaction network, we applied the gener- 
alised version of a simplified chemical network introduced 
in Ilgner & Nelson (2006a) who named this the "modi- 
fied Oppenheimer-Dalgarno model" . It is a modification of 
the reaction scheme proposed by Oppenheimer & Dalgarno 
(1974) and was obtained by adding the grain and man- 
tle chemistry. We now dropped the restriction |iT m i n | = 
l-^maxl = 2 we made in Ilgner & Nelson (2006a) and al- 
low Z m i n and Z mx to become free of choice. A full list of 
the chemical reaction schemes applied is given in Table [T] 
To keep the terminology simple, we still refer to this model 
as the " modified Oppenheimer-Dalgarno model" . 

We benefit from the fact that for our particular chem- 
ical model the numerical solution of the kinetic equations 
associated with the equilibrium state is accompanied by 
an analytical solution. The key idea needed to construct 
the analytical solution was taken from Okuzumi (2009), 
who succeeded to obtain an analytical solution for the net- 
work of Umebayshi-Nakano (1990). Our network, i.e. the 
modified Oppenheimer-Dalgarno model, and the more com- 
plex Umebayashi-Nakano network use of the same chemi- 
cal key processes. However, they differ significantly in the 
way neutral gas-phase species may stick onto grains. To fa- 
cilitate the comparison with the Umebayashi-Nakano net- 
work, we precede this section by comparing the modified 
Oppenheimer-Dalgarno with the Umebayashi-Nakano net- 
work. Both networks and the chemical models presented in 



this section are evaluated at disc midplane z/h g = un- 
der minimum mass solar nebula conditions applying Eqs. 
©, (O - (Ull), and (US]). Except for the sticking coefficient 
s c = 0.3, the ionisation rate and the elemental abundance 
of metals xm were taken from Sano et al. (2000) to facil- 
itate the comparison with previous studies. In particular, 
x~m = 7.97- 10 — 5 (5 where 8 denotes the assumed heavy metal 
gas depletion with respect to solar abundances. S = 0.02 
is regarded as the most favourable condition for retain- 
ing metals in the gas-phase. In Sec. [5] we follow a comple- 
mentary approach by making a conservative estimate with 
x M < 7.97- 10- 5 5. 

3.1. Comparison: Modified Oppenheimer-Dalgarno vs 
Umebayashi-Nakano network 

Both chemical networks are frequently applied, e.g., in 
studies on the dead zone structure of protoplanetary 
discs. In particular, under conditions typically applied 
for protoplanetary discs the ionisation degree and grain 
charging processes are closely related, cf Ilgner & Nelson 
(2006a). Tracing the dead zone structure of protoplane- 
tary discs, the chemical network associated with the mod- 
ified Oppenheimer-Dalgarno model has been compared 
with more complex chemical networks, cf Ilgner & Nelson 
(2006a) and Bai & Goodman (2009). While Ilgner & Nelson 
were concerned with conventional a-disc models, Bai and 
Goodman applied minimum mass solar nebula conditions. 
Considering a huge range of grain depletion, Bai and 
Goodman showed that the complex network generates more 
free electrons, which are less than < 2 times the values ob- 
tained for the modified Oppenheimer-Dalgarno model. 

Similar to Oppenheimer & Dalgarno (1974), the net- 
work of Umebayashi & Nakano (1990) was originally de- 
signed under dense interstellar cloud conditions. Sano et 
al. (2000) adopted the scheme of Umebayashi and Nakano 
(with some modifications) to estimate the size of the 
dead zone under minimum mass solar nebula conditions. 
Okuzumi (2009) modified the Umebayshi-Nakano network 
to account for grain charging processes on irregular porous 
dust agglomerates. Again, to keep the terminology simple, 
the chemical networks presented in Sano et al. (2000) and 
Okuzumi (2009) are referred to as the Umebayashi-Nakano 
model. 

We recall that the modified Oppenheimer-Dalgarno 
model and Umebayashi-Nakano model are quite similar in 
the way they are constructed. Grain charges up to \Z\ < 3 
are considered. The networks use the schematic compo- 
nents m, nv, M, and M+. The components M and M+ 
represent a heavy metal atom and its ionized counter- 
part. The networks differ when specifying the molecular 
ion m + . In the modified Oppenheimer-Dalgarno model, m + 
simply denotes the ion of the dominant molecule, while 
the ions 0+, Oj, OH+, 2 H+, CO+, CH+, and HCO+ con- 
tribute to m + in the scheme of the Umebayashi-Nakano 
model. Umebayashi & Nakano introduced this extension 
because molecular ions react differentially. Only a few of 
the ions do react via charge transfer, others do not. In ad- 
dition to the steady-state assumption for the neutral gas 
components CO,Mg, 2 , and O, two chemical processes 
make the chemical networks different. The reaction scheme 
of the modified Oppenheimer-Dalgarno model includes the 
thermal adsorption of neutral gas-phase particles on grain 
surfaces and its corresponding reverse desorption process, 
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Fig. 3. Equilibrium concentrations at disc midplane shown for some representative components of the Umebayashi- 
Nakano (left panel) and the modified Oppenheimer-Dalgarno network (right panel), respectively. Apart from s c = 0.3 
the parameter setup is taken from Sano et al. (2000). The elements m and M of the modified Oppenheimer-Dalgarno 
network are specified by molecular hydrogen and magnesium. The line Xi = oo at R = 10 AU is marked to aid comparison 
with the modified Oppenheimer-Dalgarno network applied to dust agglomerates, cf Figs. [4] and [5] 



the scheme of Umebayashi & Nakano does not. As Bai & 
Goodman (2009) pointed out, one can estimate the tem- 
perature range beyond the adsorption that dominates the 
desorption process, and vice versa. For heavy metal atoms 
M the corresponding temperature range is much more nar- 
row. For gas temperatures below a critical value, the grains 
are very effective at sweeping up metal atoms. This was 
demonstrated by Ilgner & Nelson (2006a) when studying 
which effect this has on the size of the dead zone in con- 
ventional a-discs. 

In a previous publication (Ilgner & Nelson 2006a) we 
have already examined the Umebayashi-Nakano network 
under minimum mass solar nebula conditions that repro- 
duces the results of Sano et al. (2000). Varying s e to 0.3, 
we have calculated the equilibrium abundances obtained for 
the modified Oppenheimer-Dalgarno and the Umebayashi- 
Nakano network, respectively. The profiles obtained for 
z/hg = are shown in the left panel (Umebayashi-Nakano 
network) and in the right panel (modified Oppenheimer- 
Dalgarno network) of Fig. [3] The components m and M 
of the modified Oppenheimer-Dalgarno model are specified 
by molecular hydrogen and magnesium with evaporation 
energies of 450 K and 5300 K, respectively. We moreover 
assumed a standard value of 1.5 x 10 15 cm~ 2 for the surface 
density of sites available on the grain surface. The compar- 
ison of the corresponding profiles revealed major changes 
in the metal ion M + abundance, which are remarkable, but 
not unexpected. According to Umebayashi & Nakano, M 
is steadily exchanging charges with the ionized molecules 
H^",C + , and HCO + via charge transfer while keeping the 
reservoir of neutral metal atoms Xoo [M] constant and equal 
to the total fractional abundance of metals^ In the modi- 
fied Oppenheimer-Dalgarno model at 135 K (R ~ 4.2 AU), 
the amount of metals adsorbed onto grain particles balances 
the gas-phase heavy metals. For lower temperatures essen- 



5 Problems related to the steady-state assumption of neu- 
tral metal atoms are discussed for the original Oppenheimer- 
Dalgarno model in Ilgner & Nelson (2006a), where it was la- 
belled model 1. 



tially no metals are left in the gas phase and hence no metal 
atoms are available to be processed to metal ions via charge 
transfer reactions. Instead, the molecular ion m + is becom- 
ing the dominant gas-phase ion for (R > 5 AU). However, 
the profiles for the ionisation fraction and the grain charges 
for the modified Oppenheimer-Dalgarno model and the 
Umebayashi-Nakano model agree very well. 

We note that metals serve as the key ingredients in 
governing the chemistry in the original Oppenheimer- 
Dalgarno model as well as in the Umebayashi-Nakano 
model|f| Comparing with the corresponding profiles ob- 
tained for the Umebayashi-Nakano network, we find that 
the amount of metal ions available is tremendously reduced 
in the modified Oppenheimer-Dalgarno network. Whether 
or not this may have significant consequences on the grain 
charging through reaction listed in Tab. Q] will be anal- 
ysed in the following subsection. 

3.2. Grain chemistry with irregular porous agglomerates 

The charging of irregular porous agglomerates was intro- 
duced by Okuzumi (2009), who applied the BCCA model 
to characterise the agglomerates. Okuzumi presented an 
analytical solution to obtain the equilibrium state of the 
ionisation degree and the excess charges Z of the dust ag- 
glomerates. The latter is well described by a Gaussian with 
a mean value of < Z > and a variance < AZ 2 > , which is 
nicely detailed in Okuzumi (2009). However, at first glance 
the analytical solution might serve to verify the numeri- 
cal solution rather than providing an a priori estimate. To 
recap: solving the corresponding equations (these are Eqs. 
(27) - (31) in Okuzumi, 2009) requires knowing the final 
composition of the ions in order to calculate the mean ion 



Under certain conditions, this may have important conse- 
quences for ongoing planet formation in discs: In the ion-electron 
plasma limit (see Subsec. fo.ip the structure of the dead zone may 
be significantly altered by the combined action of recombination 
of metal ions and turbulent transport (Ilgner & Nelson 2006b, 
2008). 
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Fig. 4. Mean charge < Z > of dust agglomerates and the standard deviation V < AZ 2 > as a function of the number 
of constituent monomers N (left panel). The upper abscissa relates ./V to the characteristic radius a c of the BCCA 
agglomerate. The electron fraction x e , the ion concentration Xi (with M + as the dominant ion) and the mean charge 
of the dust <Z> per dust agglomerate are shown at the right panel using solid, dashed, and dotted lines. The profiles 
shown are obtained at R = 10 AU and z/h g = applying the "reduced" network of the modified Oppenheimer-Dalgarno 
model with adsorption and desorption process switched off. Dbcca — 2. The lines correspond to the values obtained by 
numerical integration while the values marked by (blue) filled circles are obtained using the semi-analytical solution. 



velocity Ui and the mean gas-phase recombination rate (3. 
For a metal- rich environment, we have seen that the metal 
ion Mg + is by far the dominant ion in the Umebayashi- 
Nakano network (see previous subsection), which compen- 
sates its lower mean thermal velocity. We therefore simply 
set Hi w % g +. As already mentioned in Okuzumi (2009), 
the contributions of g(f3) are negligible for low values of the 
gas-phase recombination rate coeffcients applied. 

With respect to the electrostatic polarisation of dust 
material caused by the electric field of an approaching 
charged particle, agglomerates and spherical particles differ 
significantly. Polarisation contributes to the mean collision 
rate coefficient <av> oc J with J denoting the so-called 
"reduced rate" (Draine & Sutin, 1987): 



J 



1/2 



1-^1 



1/2 



if v = 

if ^ < (26) 



expj - %4(l+ Ut + 3v) 1/2 y else 



with 0„ s» v 1 1 (l+u^ 1 / 2 ). v — Ze/qi relates the grain charge 
Ze to the charge of the colliding particle (q = e for singly 
charged ions and q = —e for electrons) while r = akTq~ 2 
refers to the so-called " reduced temperature" . In contrast to 
spherical dust particles, agglomerates are considered to be 
hardly affected by electrostatic polarisation, see discussion 
in Okuzumi (2009). For r — > oo and \v\ 3> 1 the contribu- 
tions to J through polarisations become negligible and the 
expression for J reduces to (Draine & Sutin, 1987) 



J 



exp 



r 



if v < 
else . 



(27) 



The polarisation also effects the sticking coefficient 
s e (a, Z,T) which is defined as 



f~P e (E)a(E)eM- 



J™a(E)exp{-£}dE 



(28) 



assuming that the electron obeys a Maxwellian distribution 
at infinity We also have to specify the grain topology since 
P e depends on grain properties such as the grain surface. 
For example, Umebayashi & Nakano (1980) assumed planar 
surfaces to derive P e . We did not investigate this question 
for dust agglomerates and simply assumed that s e — 0.3 is 
independent of the grain charge. In particular, we regard 
s e as an adjustable parameter of the simulations presented, 
which agrees with Okuzumi (2009). 

When applying the formalism proposed by Okuzumi 
(2009) to the modified Oppenheimer-Dalgarno model, we 
followed a two step approach. To be able to compare our re- 
sults with Fig. [21 we applied the same parameter setup. We 
began with switching off the adsorption and desorption pro- 
cess of the modified Oppenheimer-Dalgarno network, i.e., 
^5 — kg — (referring to Tab. [T]), and ended in the ana- 
lytical solution presented in Okuzumi (2009). The results 
obtained, for example at disc midplane at R — 10 AU, are 
shown in Fig. 01 The mean value <Z> and standard devia- 
tion \/ <AZ 2 > of the excess charge Z of the dust agglomer- 
ate are shown in the left panel while the right panel reveals 
the outcome of the comparison between analytical and nu- 
merical values for the electron fraction a;[e _ ], the molecule 
and metal ion concentration x[m + ], x[M+], and the gram 
charge density approximated by <Z> Xd- The analytical 
and numerical values agree excellently, which proves that 
our implementation of Okuzumi's method is correct. 

In a final step we constructed an analytical solution that 
covers k$ks ^ 0. The rates associated with thermal adsorp- 
tion and desorption contribute exclusively to the rate equa- 
tions of the neutral gas-phase components and the mantle 
components. The kinetic equations for the non-neutral gas- 
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Fig. 5. Mean charge < Z > of dust agglomerates and the standard deviation V < AZ 2 > as a function of the number 
of constituent monomers A (left panel). The upper abscissa relates ./V to the characteristic radius a c of the BCCA 
agglomerate. The electron fraction x e , the ion concentration Xi (with M + as the dominant ion) and the mean charge 
of the dust <Z> per dust agglomerate are shown in the right panel using solid, dashed, and dotted lines. The profiles 
shown are obtained at R — 10 AU and z/h s = applying the modified Oppenheimer-Dalgarno model with adsorption and 
desorption process switched on. Dbcca — 2. The lines correspond to the values obtained by numerical integration, while 
the values marked by (blue) filled circles are obtained using the semi-analytical solution. The unlabelled (red-coloured) 
dotted lines are the corresponding profiles obtained for the "reduced" network of the modified Oppenheimer-Dalgarno 
model with adsorption and desorption process switched off, which we already showed in Fig. @] 



components remain unchanged and are (for t — > oo) 
A[m+] = k 4 N[m}/ (29) 

(k 1 N\M\ + k»Ne + v m +n 4 <a[m+ ,N d (Z)]>) 
N[M+] = k 1 N[M}N[m+}/ (30) 

(k 3 N e + v M+ n d <a[M+,N d {Z)]>) 
N e = k 4 N[m}/ (31) 

(pN t +v e -n d <a[e-,N d (Z)]>) 

with 

n d = ^N d [Z) 
z 

= k 2 N[m + ] +k 3 N[M + ], 

where A[x] denotes the number density of the component 
x while Nd(Z) is the number density of grains with charge 
Z. 

However, the (thermal) adsorption/desorption of the 
neutral gas-phase components may have an indirect effect 
on both the ionized gas-phase species and the charge state 
of the dust since this process may control the amount of 
gas-phase species available. We will demonstrate that the 
set of Eqs. (|2T)|) - (|3"Tj) and the charge balance equation 

Ni-N e + ^ZN d {Z) =0 (32) 
z 

can be solved for a single parameter <Z> if A[m] and A[M] 
serve as constants. Indeed, A[m] = const and A[M] = const 
can be extracted from the modified Oppenheimer-Dalgarno 
model. Molecular hydrogen serves as the neutral molecule 
component, which is assumed to relax towards hydrostatic 
equilibrium on a dynamical time scale t z = fl^ 1 (see sub- 
section 2.2). From the kinetic equation for the mantle com- 
ponent gM we obtained N[M] — k§Nu./{khn d + k 6 ) as- 
suming Am Ri A[M] + A[gM]; reaction kr contributes to 



<Z> > only. Solving the set of Eqs. (J29J) - ([32]), we ob- 
tained the semi-analytical solution for A^m" 1 "], A^M" 1 "], 
Aoo[e _ ], and <Z> oa (with Aoo denoting the asymptotic 
limit for t — > oo). The results obtained by numerical inte- 
gration as well as by applying the semi-analytical approach 
are presented in Fig. [5l The agreement is excellent. We 
notice that the values associated with the "reduced" net- 
work of the modified Oppenheimer-Dalgarno model (i.e., 
with adsorption and desorption switched off k§ = kg = 0) 
differ significantly from those of the "full network" with 
adsorption and desorption switched on (k 5 k 6 ^ 0). For 
A = const, the "reduced" network produces higher mean 
excess charges (associated with unlabelled red-coloured 
dotted lines in Fig. O than the same network does for 
k^ks 0. The same applies for the electron concentra- 
tion x e , < Z > Xd and the ion concentration xi with m + 
becoming the dominant ion in the modified Oppenheimer- 
Dalgarno model. Before we draw our conclusions, we inves- 
tigate this problem for the Umebayashi-Nakano network, 
too. In particular, we analysed the effect of the domi- 
nant molecular ion HCO + on the grain charging under the 
assumption that the gas-phase metal is significantly de- 
pleted through the action of the thermal adsorption. We 
therefore have modified the Umebayashi-Nakano network 
by adding the thermal adsorption/desorption process for 
MgQ The results obtained for calculations with rate co- 
efficients fcs = kg — and fcsfcg ^ 0, respectively, are 
quite different compared with the results discussed above. 
Switching on the adsorption/desorption process for met- 
als (i.e., k^k^ 0) leads to higher values for x e and Xi 
and to higher negative mean charges on dust. However, the 
changes are just in the range of a few percent and there- 
fore less severe than the changes observed in the modified 



7 Technically speaking, we included two additional rate equa- 
tions for the new components Mg and gMg. 
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Fig. 6. Mean charge < Z> of spherical dust grains and the standard deviation \/ < AZ 2 > as a function of the grain radius 
a* (left panel). The upper abscissa relates a* to a BCCA agglomerate of TV monomers assuming m d = const. The electron 
fraction x e , the ion concentration Xi (with m + as the dominant ion) and the mean charge of the dust <Z> per dust grain 
are shown in the right panel using solid, dashed, and dotted lines. The profiles shown are obtained at R = 10 AU and 
z/hg = applying the modified Oppenheimer-Dalgarno network for spherical grains. The lines correspond to the values 
obtained by numerical integration while the values marked by (blue) filled circles are obtained using the semi-analytical 
solution. To be compared with Fig.[5l which is associated with the corresponding modified Oppenheimer-Dalgarno model 
for BCCA agglomerates. 



Oppenheimer-Dalgarno network. These different features of 
the modified Oppenheimer-Dalgarno and the Umebayashi- 
Nakano network can be traced down to differences in the 
thermal velocity of dominant gas-phase ions. Under metal 
rich conditions, Mg + is the dominant ion in the modi- 
fied Oppenheimer-Dalgarno network with k$ = fc 6 = 0, 
while Hj~ for k^k^ ^ 0. Since has a much higher ther- 
mal velocity than Mg + it is quickly captured by grains. 
Applying the same conditions, HCO + becomes the domi- 
nant ion in the modified Umebayashi-Nakano network with 
a thermal velocity lower than that of Mg + , which is the 
dominant ion for k$ = k^ = 0. Because of the slightly 
reduced contributions owing to collisions between grains 
and HCO + , the values for x e and Xi are increasing. The 
drop/rise in the mean grain charge observed in the modi- 
fied Oppenheimer-Dalgarno and Umebayshi-Nakano model 
also originates in differences between the thermal veloci- 
ties of the dominant gas-phase ions: the higher the thermal 
velocity of the dominant ion, the less the negative mean 
charge on grains. Considering all this, we find indeed that 
the adsorption/desorption of metals effects the mean value 
for grain charges, the electron fraction, and the ion abun- 
dances. Whether or not the value increases/decreases de- 
pends on the ratio of the thermal velocities for the metal 
ion and the dominant molecular ion. 



3.3. Grain chemistry with spherical grains 

Instead of agglomerates we now consider spherical dust 
grains of a given radius do- Applying the modified 
Oppenheimer-Dalgarno model, we aim to determine their 
mean excess charge <Z>. If we neglect the effect of polar- 
isation, <Z> can be calculated semi-analytically following 
the same procedure detailed in the previous subsection. We 
simply repeated the procedure and obtained an identical set 



of equations to Eqs. ([23)1 - (1521) by replacing 
N d (Z) -> N dQ (Z) 

where Ndo(Z) and n d Q denote the number density of spher- 
ical grain particles with excess charge Z and the total num- 
ber density of spherical grains, respectively. The results 
obtained are shown in Fig. [6] In the left panel the mean 
grain charge and its the standard deviation are plotted as 
a function of the grain radius, while the fractional abun- 
dances obtained are shown in the right panel. We first note 
the excellent agreement between the numerical and semi- 
analytical values. We can compare the results with those 
obtained for agglomerates when identifying the grain ra- 
dius with a*, see Eq. (|18|) . Assuming monomers of size 
ao = 10~ 5 cm, a„ scales with the number N of constituent 
monomers as shown in the upper x-axis of Fig. [6] A compact 
spherical dust grain of, for example, a» = 10~ 3 cm corre- 
sponds to an agglomerate made by N = 10 6 monomers 
of size ao = 10 -5 cm. Comparing Fig. [5] with Fig. [6j we 
find that the mean charge carried by compact spheres is 
much lower than the value obtained for agglomerates. For 
N = 10 6 , it is about one order of magnitude. This is not 
caused by the changes in the projected surface area associ- 
ated with the transition between fractal and compact grain 
topology. It is because of the governing ion-electron regime 
in which the mean grain charge is proportional to the ra- 
dius, i.e. a c /a* = N 1 ^. For the same reason, the modified 
Oppenheimer-Dalgarno network applied to compact spher- 
ical grains produces a much larger amount of gas-phase ions 
m + and e _ and a much lower dust charge density <Z> x d 
than the values obtained for agglomerates. We conclude 
that dust agglomerates have a higher charge-to-mass ra- 
tio carrying more charges than the corresponding compact 
spheres. 

We completed the cycle by taking the polarisation of 
compact spherical grains into account. We cannot apply the 
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equations presented in Okuzumi (2009) though, which were 
derived under the assumption of Eq. (1271) . Substituting Eq. 
I|27p with Eq. ([2T5)) leads to different expressions for <<7di> 
and <(7de>- This may also affect the Gaussian distribution 
for Nd {Z) since the coefficient W(Z) of the correspond- 
ing first-order equation changes (cf appendix in Okuzumi, 
2009). However, we did not consider the semi-analytical ap- 
proach and performed a numerical integration of the rate 
equations taking Eq. (|2"o) into account. The values for <Z> 
and < AZ 2 > obtained under the assumption of Eq. (|2"T|) 
were used as an initial guess for the range of grain charges 
[Z m i n , Zmax] considered. If necessary, the charge range was 
modified and the integration was repeated again until the 
calculated mean value <Z> fitted the peak value of the 
Gaussian profile for N<io{Z). The values obtained by nu- 
merical integration match the analytical values associated 
with the non-polarised case very well. The corresponding 
profiles show no distinguishing features when applying the 
logarithmic scaling of Fig. |51 Changing to linear scales, we 
observe for sizes a < 10~ 4 cm that the analytical values (i.e. 
the model neglecting the polarisation) obtained for < Z > 
and < AZ 2 y 1 / 2 are slightly lower than those obtained by 
numerical integration (model with polarisation). We con- 
clude that at least for our particular disc model polarisation 
has only a minor effect on grain charging. 

4. lonisation rates 

Because of the simplicity of the chemical network applied, 
all processes that might result in ionising disc material 
could not be treated individually. Instead, the contributions 
of the dominant processes were combined into an effective 
ionisation rate £ with &4 = £ (see Tab.Q]). 

We included the ionisation of the dominant molecule 
by incident X-rays that originate in the corona of the cen- 
tral T Tauri star. The corresponding X-ray ionisation rates 
Cxray can be obtained from radiative transfer calculations. 
The results of these calculations are presented in Igea & 
Glassgold (1999), who modelled the transport of stellar, 
low-energy X-rays (< 20 keV). The ionisation rates are 
published primarily for the minimum solar nebula model of 
Hayashi et al. (1985), which corresponds to the static disc 
model discussed above. Taking the photoelectric absorption 
and Compton scattering into account, Igea & Glassgold 
specified absorption-dominated and scattering-dominated 
regimes, respectively. Bai & Goodman (2009) presented a 
best fit to the X-ray ionisation rates of Igea & Glassgold 
(1999), which facilitates employing the results of Igea & 
Glassgold (1999). 

Igea & Glassgold (1999) presented ionisation rates ob- 
tained for power indices p s of the gas surface density E g 
except for p s — —3/2. For a specified orbital radius at 
R = 1 AU they report that the deviation in the X-ray 
ionisation rate is generally less than a factor of two. Igea 
& Glassgold stated that the ionisation rates plotted vs the 
vertical column density "manifest a universal form" that 
is independent of the surface density profiles considered. 
Glassgold et al. (2000) note that this "scaling is also found 
at other disk radii and is insensitive to changes in disk pa- 
rameters such as disk mass, surface density, and slope q."0 
However, X-ray ionisation rates obtained from X-ray trans- 

8 In Glassgold et al. (2000) q refers to the power index asso- 
ciated with the radial variation of the gas density. 



fer calculations have not been made publicly available yet 
for profiles of the gas surface density different from the 
static model of Hayashi et al. (1985). 

Instead of a complex X-ray radiative transfer model, we 
attempted to apply a simple ray tracing model neglecting 
scattering effects. The techniques employed are described in 
detail in, e.g., Fromang et al. (2002) and Ilgner & Nelson 
(2006a). For static and stationary disc models of minimum 
mass solar nebula type, these calculations can be performed 
in straight line because of the non-discrete nature of the 
local gas density g g . We find significant differences when 
comparing the X-ray ionisation rates we obtained with the 
rates of Igea & Glassgold (1999). Similar findings are re- 
ported in Bai & Goodman (2009). The effect of scatter- 
ing may explain these differences, but additional radiative 
transfer calculations would make valuable contributions. 

Acknowledging the potential X-ray scattering may have 
on the ionisation rate, we applied X-ray ionisation rates 
based on the simulations of Igea & Glassgold (1999). In par- 
ticular, we applied the fitting formula associated with Eq. 
(21) of Bai & Goodman (2009), which is valid for plasma 
temperature k^T = 3 keV. However, because of the prob- 
lem discussed above we regard the X-ray ionisation rates 
Cxray as an order of magnitude estimate. 

We also considered cosmic ray particles as a source for 
ionisation with ionisation rates £cr given by eq. (19) of 
Sano et al. (2000). We adopted the nominal value Co = 
10 -17 s _1 of galactic cosmic rays associated with unatten- 
uated penetration of cosmic ray particles. We neglected the 
minor contributions from the decay of radioactive elements. 
Taking the X-ray and the cosmic ray ionisation rates at 
face value, X-rays dominate the cosmic ray ionisation for 
regions close to the disc surface z/h g = 3. At higher op- 
tical depths, ionisation through cosmic rays is the domi- 
nant source of ionisation. We regard X-rays as our stan- 
dard ionisation source while ionisation through cosmic rays 
is restricted to specific disc regions. For inner disc regions 
R < R\ we excluded cosmic rays and set £ = Cxray We fur- 
thermore assumed that the shielding effect of the T Tauri 
winds peters out between Ri < R < R 2 . We assumed a 
simple exponential decay 

to mimic this effect for that particular transition region 
from X-ray dominated to cosmic ray dominated regions 
with Ci = CxrayCRi, z) and ( 2 = Ccr(-R2, z). 

5. Numerical method 

Although we applied a two-dimensional geometry (R,z), 
the dynamics of the gas-dust disc refers to a one- 
dimensional description. Therefore, no serious attempt was 
made to couple the dynamical and the chemical evolution of 
the gas-dust disc by means of conventional operator split- 
ting techniques. The coupling is actually not required since 
the grain charging operates on a much shorter time scale 
than the disc dynamics. However, we aimed to simulate 
the disc dynamics to halt the simulation at different evolu- 
tionary stages. Assuming that the local disc structure is set 
up instantaneously, the kinetic equations are integrated nu- 
merically for t = 10 4 yr at each point of the two-dimensional 
domain. In parallel we determined the limiting behaviour 
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(i.e., for t — > oo) following the semi-analytical approach. 

In order to simulate the disc dynamics we continued by 
applying the so-called "method of line" (MOL) approach 
as we did, e.g., in Ilgner & Nelson (2006b), to solve the set 
of coupled advection-diffusion equations (U]) and (|5|). We 
applied a radial grid {Ri : i = 1, • • • , n^} with an equidis- 
tant mesh size near the inner boundary and non-equidistant 
mesh sizes elsewhere. The radial domain R = [0.1, 10 3 ] AU 
was discretised using 71r = 222 grid cells. The transport 
terms were discretised in the flux form to secure the conser- 
vation of mass. We applied the concept of staggered mesh 
on which the scalars and vectors representing the discrete 
values of the independent variables are centred at different 
locations. Scalars are located at zone centers while vectors 
are defined at zone edges. We limited the maximum abso- 
lute step size corresponding to a Courant-Friedrichs-Lewy 
(CFL) number of 0.5. 

We applied the Crank-Nicholsen scheme when discretis- 
ing the diffusion operator in Eq. ((SJ, as we did in Ilgner 
& Nelson (2006b). Regarding the advection operator, we 
employed a nonlinear advection discretisation scheme that 
prevents negative values for the individual mass densities. 
In particular, we employed an upwind biased discretisation 
scheme with flux limiting. We opted for the flux limiter 
function proposed by Koren (1993), which - depending on 
the smoothness of the profile for the quantity to be advected 
- is associated with the accuracy of a third-order scheme. 

The kinetic equations associated with the chemical 
model are given by a set of ordinary differential equations 
(ODE). We employed stiff ODE integrators to obtain sta- 
ble numerical solutions. Stiffness was introduced for both 
physical and numerical reasons because of, e.g., the huge 
range of time scales for the chemistry and because of the 
semi-discretisation associated with the MOL approach. In 
particular, we used standard stiff ODE solvers based on 
linear multistep methods, namely the backward differenti- 
ation formulas. 

Concerning the semi-analytic approach, we employed 
the procedure proposed by Muller (1956) to find the roots 
of the polynomial associated with Eqs. (l29l) - ([32]). 



6. Results 

In this section we present the results of simulations that 
examined the modified Oppcnhcimcr-Dalgarno network un- 
der non-static disc conditions. We examined the modified 
Oppenheimer-Dalgarno network under stationary disc con- 
ditions which we later used for a comparison with the 
results obtained for non-stationary discs. We evolved the 
disc chemistry for t — 10 4 yr. However, by applying the 
semi-analytic approach discussed in Sec. 02 we are in the 
favourable position to calculate the long-term behaviour of 
the chemical models a priori. That is why we can prove the 
time interval needed to establish an equilibrium solution, 
which is (under the conditions applied) shorter than the 
dynamical time scale £k- 

Apart from the radius a* and a c , respectively, we kept 
all other parameters fixed: M* = M@, [i — 2.33, a — 10 -3 , 
g p = 1.0 gem" 3 , Sc = 1, Vg = 3.33 x 10~ 2 AU, and \M d \ = 
10" 10 Moyr" 1 . We adopted values L x = 10 29 erg s _1 
and fceTx = 3 keV for the total X-ray luminosity and 
the plasma temperature while the transition between X- 
ray dominated and cosmic ray dominated ionisation region 



was fixed to R 1 = 25 AU and R 2 = 30 AU (see Sec. g|). 
The values for i?i and R 2 were motivated by Perez-Becker 
& Chiang (2011), who discussed the problem of so-called 
" sideways cosmic rays" . Regarding the parameter setup as- 
sociated with the chemical model we fixed the sticking co- 
efficients s = 0.3 and Si — 1.0, which is in line with the val- 
ues Okuzumi (2009) applied. As we did in Ilgner & Nelson 
(2006a), we assumed the standard value of 1.5 x 10 15 cm 2 for 
the surface density of sites, while the value of energy AE$ 
transferred to the grain particle caused by lattice vibration 
was approximated by 2.0 x 10~ 3 eV. We approximated the 
dissociation energy D for neutral gas-phase components by 
its binding energy for physical adsorption. Complementary 
to the results discussed in Sec. 131 we applied a conservative 
estimate of % = 10 -11 in this section. 

We are aware that the disc model we applied serves 
as a toy model and does not cover the complexity of the 
dust-gas dynamics, a = 10~ 3 is indeed a very naive ap- 
proximation for the turbulence structure in discs and, in 
particular, in dead zones. However, recent simulations have 
shown that sound waves propagating into dead zones pro- 
mote the diffusion of dust particles (Okuzumi & Hirose 
2011). Lacking reliable models, we excluded grain size 
distributions. Instead we considered mono-sized particles 
and agglomerates and set limits on a c and a*. We chose 
a* = [4.64- 10~ 5 , 10~ 3 ] cm, which corresponds to agglomer- 
ates made of N = [10 2 , 10 6 ] monomers of size ao = 10~ 5 cm. 
This limitation is motivated by recent studies of Okuzumi 
et al. (201 lab). Evaluating a set of growth criteria for 
BCCA agglomerates, they made predictions on the location 
of so-called " frozen zones" . In frozen zones the electro-static 
barrier inhibits the agglomerates to growth further. The 
authors showed that planet forming regions are associated 
with frozen zones. The inferred sizes of the agglomerates 
are [4 • 10~ 5 , 3 • 10" 2 ] cm, depending on the local position. 

6.1. Stationary disc 

The assumption of stationary discs, i.e., if the mass flow is 
constant with the orbital radius, corresponds to a gas sur- 
face density profile given in Eq. © with Eo = 350 gcm~ 2 . 
One primary aim of this work is to determine to what ex- 
tend compact spherical grains and dust agglomerates may 
allow for different grain charges Z. We already answered 
this question for a fixed dust-to-gas ratio Ed/S g = 10~ 2 
(see the preceding section [3]); in this section we conduct a 
similar exercise for a fixed mass flow of dust material 

We already know that the population N d (Z) (and 
Nao(Z)) follows a Gaussian and therefore is completely 
characterised by the mean value < Z > and the variance 
<AZ 2 >. In particular we compare compact spherical grains 
and dust agglomerates of the same mass, which constrains 
the radius a* of the spherical grain. For BCCA agglomer- 
ates (with -Dbcca = 2) relation Eq. (|18D holds. 
Figure shows the result of the one-to-one comparison. 
In the left panel of Fig. [7] the mean value <Z> and the 
standard deviation V < AZ 2 > obtained for dust agglom- 
erates at different altitudes are plotted against the orbital 
radius R. The agglomerates consist of N monomers of size 
ao = 10~ 5 cm; the size of the agglomerate varies from 
N = 10 6 (top) to N = 10 2 (bottom ) passing N = 10 4 . 
The profiles of <Z> and V < AZ 2 > obtained for the cor- 
responding compact spherical grains are shown in the right 
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Fig. 7. Mean charge < Z > and the standard deviation V < AZ 2 > associated with two different dust topologies as a 
function of the orbital radius R with |Md| = 10 -10 M* yr -1 . The profiles shown in the left panel refer to BCCA 
agglomerates with N = 10 6 (top), N = 10 4 (middle), and N = 10 2 (bottom) monomers of ao = 10~ 5 cm. The results 
obtained for the corresponding compact spheres are shown in the right panel. Here, the radius a* of the compact sphere 
varies from a* = 1CP 3 cm (top) to a* = 4.64 • 10~ 5 cm (bottom) passing a* = 2.15 • 10~ 4 cm. The lines correspond 
to the values obtained by numerical integration while the values marked by (blue) filled circles are obtained using the 
semi-analytical solution. 



panel with a* — 10 3 cm (top), a* = 2.15-10 4 cm (mid- 
dle), and a* = 4.64 • 10" 5 cm (bottom). While the solid 
and dotted lines refer to the values obtained by numerical 
integration, the values marked with (blue) filled circles cor- 
respond to the semi-analytical solution. 

The first thing we comment on is the significant de- 
viation of the analytical value of v < AZ 2 > from the nu- 



merical value at disc midplane z/h s = for R < 1 AU. 
Each numerical value shown in Fig. [7] corresponds to the 
numerical solution of the ODE system for At = 10 4 yr. 
We confirmed by numerical integration that < AZ 2 > ap- 
proaches its equilibrium state beyond At = 10 4 yr because 
of the low ionisation rates in this particular region. For ag- 
glomerates with N = 10 6 we obtained, for example, at disc 
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Fig. 8. Population Nd(Z) and Ndo(Z) (measured in parti- 
cle concentration), respectively, vs charge excess Z eval- 
uated at disc midplane at R = 10 AU. The Gaussian 
distribution of Nd(Z)/n s with a peak around Z = —8 
refers to agglomerates made of N = 10 6 monomers of size 
ao = 1CP 5 cm while the profile with a peak at Z = —70 is 
obtained for compact spheres with a* = 10~ 3 cm. The dot- 
ted and the dashes lines mark <Z> and <Z> ±\/ <AZ 2 >, 
respectively. The solid lines correspond to the values ob- 
tained by numerical integration while the values marked by 
(blue) filled circles are obtained using the semi-analytical 
solution. 



midplane at R = 0.7 AU \/<AZ 2 > oc = 31.54 instead of 
y/ < AZ 2 > = 18.97. However, we also remark that because 
of <Z> —> 0~ with decreasing R, the semi-analytical ap- 
proach may fail to comply with the assumption Z < 0, 
which was made to construct the semi-analytical solution. 

We found for the population of BCCA agglomerates 
that the grains follow the same population N&{Z) for re- 
gions R > 30 AU independent of height z/h g . That is partly 
because cosmic ray particles become the dominant source 
of ionisation beyond R = 25 AU and therefore can pene- 
trate easily towards disc midplane because of the low gas 
surface densities. However, we stress the more important 
matter of the ion-electron regime associated with these re- 
gions. Here Nd(Z) is determined by the grain size and the 
gas temperature, which is assumed to be independent of z. 

Because T$ oc ao, the three simulations presented in the 
left panel of Fig. [7] correspond to the same radial profile 
of Ed- We find that the mean value and the standard de- 
viation decrease by reducing N and < Z > oc A^ for the 
ion-dust regime in particular (Okuzumi 2009). The profiles 
in the right panel of Fig. [7] show similar features. 

In a next step we examined the population of BCCA ag- 
glomerates and compact spherical grains with respect to 
their excess charge. To recap: The populations Nd(Z) and 
Ndo(Z) are assumed to obey a Gaussian 



N d (Z) 



n d f l (Z-<Z>) 

V2ir <AZ 2 > GXP i 2 <AZ 2 > 



(34) 



Recalling the substantial changes in <Z> and V < AZ 2 > 
(cf Fig. [7]), we expect significant changes with respect to 
their population. This is indeed what we observe in Fig. [5] 
There we plot the population Nd(Z) and Ndo(Z) (in terms 



of particle concentration) at disc midplane at R = 10 AU 
for BCCA agglomerates with N = 10 6 and compact spheres 
of a* = 10 -3 cm. Again, the solid lines refer to the val- 
ues obtained by numerical integration while the values 
marked with (blue) filled circles correspond to the semi- 
analytical solution. Comparing with the corresponding pro- 
file for BCCA agglomerate, we find that for a given Md 
i) the mean value obtained for spherical grains is shifted 
to lower values while ii) the distribution becomes much 
narrower. Differences regarding max[xd(Z)] are caused by 
small changes in Ed because of slightly different values for 
T s . 

We compared the results obtained for BCCA agglomer- 
ates and the corresponding compact spheres with respect 
to the fractional abundances of the charge carrier Xi, x e , 
and <Z> Xd- The fractional abundances are shown in the 
left and right panels of Fig. U3 which correspond to those of 
Fig. [7] The values obtained by numerical integration and 
the semi-analytical solution (marked with (blue) filled cir- 
cles) agree excellently. An important discriminant is given 
by xi ss x c separating the ion-electron plasma limit from 
the ion-dust plasma limit. The ion-dust plasma limit corre- 
sponds to X[ ^> x c . The results obtained for BCCA agglom- 
erates indicate that the planet forming regions are always 
associated with the ion-dust regime. The profiles for x c , a;;, 
and < Z > Xd at disc midplane are also interesting. The 
profiles for N = 10 2 exactly fit those obtained for N = 10 4 
and Af = 10 6 , and so does N = 10 4 and N = 10 6 . The 
reason for this is simple since the cumulative projected sur- 
face area of all BCCA agglomerates is independent of N. 
The situation changes slightly when comparing the pro- 
files obtained for compact spheres, as shown in the right 
panel of Fig. [9] Here, the transition from the ion-dust to 
the ion-electron regime becomes apparent in the profiles at 
z/h g = 2 in particular. The corresponding profiles at disc 
midplane also change systematically with increasing a* . 



6.2. Non-stationary disc 

We already analysed the grain charging i) under static 
disc conditions (see section [3]) and ii) under stationary disc 
conditions, see the preceding subsection. Comparing these 
with our results obtained under stationary disc conditions, 
we now explore the impact of the gas-dust dynamics on 
the grain charging. We recall the gas-dust dynamics that 
we analysed in section [2] Compact spherical grains with 
a* < 10~ 2 cm are tightly coupled to the gas because of 
Ts <C 1. The turbulent mixing therefore dominates the ad- 
vection of grains towards the central object. For dust parti- 
cles treated as compact spheres we therefore expect that the 
evolving gas-dust disc dynamics causes minor changes in 
<Z> 7 <V AZ 2 >, and the fractional abundances compared 
with the corresponding profiles obtained under stationary 
disc conditions. Regarding the BCCA agglomerates, we al- 
ready know that the agglomerates with different N obey 
the same gas-dust disc dynamics. In particular, we find 
that the BCCA agglomerates are well coupled with the gas 
with no substantial inward migration (cf. right panel of Fig. 
H}. We conclude that the BCCA agglomerates considered 
(with A" < 10 6 ) and the corresponding compact spheres 
(with a» < 10~ 3 cm) basically follow the similar dynamical 
evolution. 

We repeated the simulations of section [5] and let the 
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Fig. 9. Fractional abundances of the charge carrier associated with two different dust topologies as a function of the 
orbital radius R with |M<i| = 10 -10 M* yr -1 . The profiles shown in the left panel refer to BCCA agglomerates with 
N = 10 6 (top), N = 10 4 (middle), and N = 10 2 (bottom) monomers of o,q = cm. The results obtained for the 

corresponding compact spheres are shown in the right panel. Here, the radius a* of the compact sphere varies from 
a* = 10~ 3 cm (top) to a* = 4.64 • 10 -5 cm (bottom) passing a* = 2.15 • 10~ 4 cm. The lines correspond to the values 
obtained by numerical integration while the values marked by (blue) filled circles are obtained using the semi-analytical 
solution, see the discussion in the text. The (red-coloured) dashed lines refer to <Z> 



gas-dust dynamics evolve. After t = 10 6 yr we stopped the 
calculation and examined the grain charging under local 
disc conditions. The profiles for <Z> are shown in Fig. [TU] 
assuming BCCA agglomerates (left panel) and the corre- 
sponding compact spheres (right panel). In particular, we 
focused on grain charging of agglomerates made of N = 10 2 
(bottom) and N = 10 6 (top) constituent monomers and the 



corresponding compact spheres of a» = 10 -3 cm (top) and 
a* = 4.64 ■ 1CP 5 cm (bottom). We also compared the pro- 
files obtained after t — 10 6 yr of dynamical evolution with 
the initial profiles at t — yr (marked with dotted lines in 
Fig. ITU|). The results confirm our expectations demonstrat- 
ing that the disc dynamics causes minor changes in the 
grain charging. However, we are aware that the situation 
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Fig. 10. Snapshots at t = yr and t = 10 6 yr of dynamical disc evolution showing the mean charge <Z> associated 
with two different dust topologies as a function of the orbital radius R. The profiles shown in the left panel refer to 
BCCA agglomerates with N = 10 6 (top) and N = 10 2 (bottom) monomers of ao = 1CP 5 cm. The results obtained for 
the corresponding compact spheres are shown in the right panel. Here, the radius a* of the compact sphere varies from 
a* = 1CP 3 cm (top) to a* = 4.64 • 1CP 5 cm (bottom). The dashed lines are associated with the profiles evaluated after 
t = 10 6 yr disc evolution; to aid comparison with the initial profiles we superimposed the corresponding profiles of Fig. 
[7] using dotted lines. As we did in the preceding figures, we mark the semi-analytical solution with (blue) filled circles. 



may change if a more appropriate X-ray ionisation rate is 
taken, see discussion in section [4] 

7. Summary 

We have presented calculations of the grain charging under 
conditions that mimic different stages of protoplanetary 
disc evolution. Instead of parametrising the dust-to-gas 
ratio, we inferred the value for Sd(i?)/S g (i?) from the 
evolution of the gas-dust disc. In particular, we applied 
the disc model of Takeuchi & Lin (2002, 2005) and 
content ourself with order-of-magnitude estimates. We 
considered collisional charging as the dominant process to 
determine the charge state of grains. For that purpose, 
we generalised the modified Oppenheimer-Dalgarno model 
introduced in Ilgner & Nelson (2006a) to account for higher 
grain charges. Based on that simple chemical network, 
we examined the grain charging for two different types 
of grain topology: compact spherical grains and fractal 
agglomerates of Df — 2. Our main conclusions are: 



1. The effect that thermal adsorption/desorption of metals 
has on grain charging depends on the mass of the domi- 
nant molecular gas-phase ion. If its mass is heavier than 
that of metal, the inclusion of the thermal adsorption of 
metals results in high negative excess charges on grains 
on average. Less negative charge excess on average is 
observed for molecular ions lighter than metals. 

2. We extended the semi-analytical method proposed 
by Okuzumi (2009), which allowed us to determine 
steady-state solutions for the mean grain charge, elec- 
tron and ion abundances associated with the modified 
Oppenheimer-Dalgarno model. The semi-analytical so- 
lutions were derived for both grain topologies: compact 
grains and fractal aggregates. 

3. The results obtained confirm that dust agglomerates 
have a higher charge-to-mass ratio than the correspond- 
ing compact spheres. 

4. We found that reducing the number N of constituent 
monomers causes a drop in the mean value < Z > of 
the grain charge and the standard deviation <V AZ 2 >. 
Under conditions valid for the ion-dust plasma limit 
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(i.e., Xi 3> x D ) the profiles for <Z> x<x, x c , and x; 
remain unchanged with varying N. This is because the 
cumulative projected surface area of all BCCA agglom- 
erates is independent of N. 

5. The results obtained by switching from one type of grain 
topology (fractal agglomerates) to another (compact 
spheres) reveal that for compact spherical grains <Z> 
is shifted towards more negative values while <V AZ 2 > 
decreases. 

6. The results obtained for agglomerate sizes a c = 
[1CP 4 , 10~ 2 ] cm indicate that the grain charging of 
BCCA agglomerates is governed by the ion-dust plasma 
limit. Transitions from the ion-dust to the ion-electron 
plasma limit are observed for compact spheres depend- 
ing on altitude z/h g . 

7. Another conclusion is that grain charging in a non- 
stationary disc environment is expected to lead to sim- 
ilar results as long as the effective ionisation rate and 
the temperature of the gas match the stationary values. 

We regard the drastically simplified description of the tem- 
peratures of the gas and the dust mentioned above as a po- 
tentially serious omission of our model. Working out more 
realistic conditions may result in significant changes in the 
local disc structure and might be a potential problem for 
future investigations. 
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Appendix A: Initial-boundary-value problems 
concerning the gas-dust dynamics 

We adopted the disc model of Takeuchi & Lin (2005) to sim- 
ulate the dynamics of the gas-dust disc. The corresponding 
set of coupled partial differential equations ^ and §5§ are 
subjected to initial and boundary conditions. We already 
know that for different sets of the dust parameters (a, g p ) 
the dynamics of the dust disc remains unchanged if the 
parameter sets correspond to the same stopping time Ts- 
Takeuchi & Lin (2005) have introduced a non-dimensional 
parameter A to retain the same dust evolution for different 
values of the initial gas mass. In the appendix we study 
the effect of the boundary conditions and the initial profile 
for Ed on the dynamics of the gas-dust disc adopting the 
fiducial model of Takeuchi & Lin (2005). 

We therefore solved Eqs. ([U and (0 for dust sizes 
of do = 1 mm assuming E g oc R^ 1 for t = to while 
Sd/Sg = 10~ 2 is truncated at R = 100 AU. Concerning the 
boundary conditions for E g and Ed, Takeuchi and Lin ap- 
plied zero torque boundary conditions at i?i nnC r = 5/6 AU 
and i? utcr = 10 3 AU. The profiles we obtained are shown 
in the left panel of Fig. IA.11 which should be compared 
with the upper panel of figure 4 of Takeuchi & Lin (2005). 
They correspond very well. However, we briefly recall the 
initial profile for the gas surface density, which is marked 
in Fig. IA.1I using (red-coloured) dotted lines. It represents 
the steady-state solution E g cx R~ x discussed in Sec. 12.11 
which was obtained under the assumption v t cx R. Applying 
the same power law v t oc R, Lynden-Bell & Pringle (1974) 
presented a similarity solution of Eq. ^ with dimension- 
less variables X = R/Rq and r = t/to + 1. In the limit of 
X<1 the expression reduces to E g cx X -1 . At this point, 
the temporal change of the gas surface density shown in 
the left panel of Fig. IA.1I becomes important. In partic- 
ular, the positive gradients in the gas surface density at 
R < 2 AU clearly indicate an inner boundary problem as- 
sociated with the zero torque boundary conditions applied 
at R = 5/6 AU. 

In order to mimic the dynamics of the gas-dust disc at 
the very inner disc regions more precisely, we set the inner 
boundary at i?i nn0 r = 0.1 AU. We also replaced the zero 
torque boundary conditions with predefined analytical val- 
ues for E g (i) applying the self-similarity solution given in 
Eq. (|15p . The most appropriate value for the constant a 
was applied to approximate E g (i) in the outer disc regions. 
Applying the improved boundary conditions, we then re- 
peated the simulation and obtained the profiles shown in 
the right panel of Fig. IA.1I Comparing with the correspond- 
ing profiles in the left panel, we find that the time evolution 
of E g at the outer disc regions remains unchanged while we 
observe E g oc R^ 1 in the inner disc regions. The changes 
in Ed at late evolutionary stages are remarkable. While in 
Takeuchi & Lin (2005) the entire dust disc has vanished 
for t = 10 6 yr, the improved boundary conditions cause a 
significant slow-down the rapid accretion of mm-sized com- 
pact spherical grain particles onto the star. 

We proceeded by reassessing the assumption Ed/E g = 
const at t = to that Takeuchi & Lin (2005) applied. We 
already mentioned the steady-state solutions © and ©. 
For small dust sizes a < 10~ 3 cm, Ed/E g = const corre- 
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Fig. A.l. Radial profile of the surface densities of the gas and the dust particles at different time snaps t = 10 5 ,3 ■ 10 5 
and 10 6 yr with ao — 1 mm and g p — 1.0 gcm~ 3 . The parameter are a = 10 -3 and £d/E g = 10 -2 at t/tpc = 0. The solid 
lines correspond to the dust surface density while the gas surface density is shown using the dashed lines. The dotted 
line refers to £ g ~ R^ 1 at t/t v = 0. The profiles shown in the left panel are obtained using the boundary conditions of 
Takeuchi & Lin (2005). Moving the inner boundary to smaller R and applying an analytical prescription for £ g at the 
inner and outer boundary cause the profiles to change as shown in the right panel. 
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Fig. A. 2. Radial profile of the surface densities of the gas and the dust at different time snaps t — 10 5 , 3 • 10 5 and 10 6 yr. 
The solid lines correspond to the dust surface density while the gas surface density is shown using the dashed lines. 
The (red-coloured) dotted line refers to S g ~ R~ x and E d = M d /(2irR < v R}d >) at t/t K = 0. \M d \ = lO" lo M /yr, 
a = 10~ 3 , and g p — 1.0 gem -3 . The profiles shown in the left panel are obtained under the assumption of spherical 
grains with ao = 1 mm. The corresponding profiles obtained for agglomerates with -Dbcca = 2 and N = 10 12 monomers 
of ao = 10~ 5 cm are shown in the right panel. 



sponds very well to © and for larger grains it does not. 
In stationary discs larger particles migrate inwards much 
faster, causing a size fractionation reported in Takeuchi & 
Lin (2002). Hence we repeated the calculation replacing the 
initial assumption Sd/S g = const by \Md\ = const. For the 
sake of convenience we applied the value of Takeuchi & Lin 
(2002) M d = lO _lo M yr _1 , which is one order of mag- 
nitude lower than the corresponding mass flux of the gas 
Af g = 37w£ g - 2.6 x 10~ 9 M Q yr^ 1 . The results obtained 
are shown in the left panel of Fig. IA.2I The initial profiles 
of E d and S g are shown using (red-coloured) dotted lines 
and confirm the size fractionation and Sd/S g ^ const at 
t = t . 



The grain particles are now regarded as BCCA agglom- 
erates with a characteristic radius r c defined by Eq. (|16p . 
We can relate the single agglomerate made of N monomers 
of size ao to a compact sphere of the same mass apply- 
ing Eq. (IT%1) . According to Eq. (Tig)) , a compact sphere of 
a* = 10 _1 cm corresponds to an agglomerate of N = 10 12 
monomers with ao = 10~ 5 cm. The surface density profiles 
obtained for such a dust agglomerate are shown in the right 
panel of Fig. IA.2I We observe that the dust is tightly cou- 
pled to the gas dynamics expanding because of turbulent 
mixing towards larger orbital radii R. 
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